Created by Ronnie Dutta and Ed Keys
Email: edward.keys15@imperial.ac.uk
HTML Version (This will be a link)
In [ ]:
# Import various packages and modules to be used
%matplotlib notebook
import numpy as np
import cmath
import matplotlib.pyplot as plt
import matplotlib.animation as animation
For $\theta_i > \theta_{crit}$ we can write
$$\cos{\theta_t} = (1 - \sin^2{\theta_t})^{1/2} = (1 - (\eta_1 / \eta_2)^2 \sin^2{\theta_i})^{1/2} = ib$$Then
$$r_s = \frac{(\eta_1 / \eta_2) \cos{\theta_i} - \cos{\theta_t}}{(\eta_1 / \eta_2) \cos{\theta_i} + \cos{\theta_t}} = \frac{a_s - ib}{a_s + ib}$$where $a_s = (\eta_1 / \eta_2) \cos{\theta_i}$. Similarly,
$$r_p = -\frac{(\eta_2 / \eta_1) \cos{\theta_i} - \cos{\theta_t}}{(\eta_2 / \eta_1) \cos{\theta_i} + \cos{\theta_t}} = -\frac{a_p - ib}{a_p + ib}$$where $a_p = (\eta_2 / \eta_1) \cos{\theta_i}$.
Both have $\left|r\right| = \frac{\left|a - ib\right|}{\left|a + ib\right|} = 1$
The transmitted wave is
$$E_t = E_{t0} e^{i (\mathbf{k}_t \cdot \mathbf{r} - \omega t)} = E_{t0} e^{i (k_t \sin{\theta_t} x + k_t \cos{\theta_t} y - \omega t)}$$Using $k_t = \frac{\eta_2 \omega}{c}$,
$$k_t \sin{\theta_t} = \frac{\eta_2 \omega}{c} \frac{\eta_1}{\eta_2} \sin{\theta_i} = \frac{\eta_1 \omega}{c} \sin{\theta_i} = k_i \sin{\theta_i}$$$$\Rightarrow E_t = E_{t0} e^{i \left(k_i \sin{\theta_i} x + i k_t b y - \omega t\right)} = E_{t0} e^{i \left(k_i \sin{\theta_i} x - \omega t\right)} e^{-k_t b y}$$This solution is a surface wave that decays away from the interface with a scale distance of $\frac{1}{k_t b} = \frac{c}{\eta_2 \omega b}$
Python Now we would like to be able to visualise what's going on so let's write some code to be able to see what's happening with a contour plot. This code also relies on knowledge of Classe so if you're not familiar with them have a read through this.
In [ ]:
class Boundary:
def __init__(self, angle, n1, n2, polarisation, interference=True, frustration=False):
"""
Creates a list of frames that progress in time
Args:
angle (float) - [degs] {0 to 90}
# out (bool) - True if outgoing, False if incoming (to the origin)
# E_0 (float) - Magnitude of the E field
polarisation (str) - {'s' or 'p'} Whether E or B is parallel to the boundary
# w (float) - [rad s^-1] Angular frequency (default: green light)
n1 (float) - The incident material's refractive index
n2 (float) - The second material's refractive index
interference
"""
self.theta = np.deg2rad(angle)
self.n1 = n1
self.n2 = n2
self.theta_crit = self.find_critical_angle()
if polarisation == "s" or polarisation == "p":
self.polarisation = polarisation
else:
raise Exception('Polarisation argument must be "s" or "p"')
self.frame = 0
self.num_frames = 100
if self.theta > self.theta_crit:
self.frustration = frustration
else:
self.frustration = False
self.graph_dim = 100
self.incident = self.create_wave(theta=self.theta, amplitude=1, material=1)
self.transmitted = self.transmit()
self.frustrated = self.create_wave(theta=self.theta, amplitude=1, material=3)
if interference and self.n1 != self.n2:
self.reflected = self.reflect()
self.incident["zz"] += self.reflected["zz"]
def create_wave(self, theta, amplitude, material, reverse_phase=False):
"""
Args:
theta (float) - [rads] {0 to π/2}
amplitude
material (int) - {1 or 2} Whether the wave is in material 1 or 2
reverse_phase (bool)
"""
if material == 1:
x = np.linspace(-1, 0, self.graph_dim/2)
n = self.n1
elif material == 2:
if self.frustration:
x = np.linspace(0, 0.1, self.graph_dim/20)
n = self.n2
else:
x = np.linspace(0, 1, self.graph_dim/2)
n = self.n2
elif material == 3:
x = np.linspace(0.1, 1.0, 9*self.graph_dim/20 )
n = self.n1
else:
raise ValueError("arg 'material' must be 1, 2 or 3")
y = np.linspace(-1, 1, self.graph_dim)
xx, yy = np.meshgrid(x, y)
k_x = np.cos(theta) * n
k_y = -np.sin(theta) * n
decay = 1
shift = 0
if material == 3:
shift = 0.1
k_t = 16*np.pi*self.n2 # this is from the wavelength being 1/8
cos_theta_t = self.cos_snell(self.n1, self.n2, self.theta)
b = -1j*cos_theta_t
decay = 1/(k_t * b)
decay = decay.real
print('decay is ', decay)
if reverse_phase == False:
zz = np.array([decay*amplitude * np.exp(8j*np.pi * (k_x*(xx-shift) + k_y*yy - phase)).real for phase in np.linspace(0, 0.25, self.num_frames)])
elif reverse_phase == True:
zz = np.array([decay*amplitude * np.exp(8j*np.pi * (k_x*(xx-shift) + k_y*yy + phase)).real for phase in np.linspace(0, 0.25, self.num_frames)])
else:
raise TypeError("arg `reverse_phase` must be of type('bool')")
return {"x": x, "y": y, "zz": zz}
def transmit(self):
cos_theta_i = np.cos(self.theta)
cos_theta_t = self.cos_snell(self.n1, self.n2, self.theta)
plot_theta_t = cmath.acos(cos_theta_t)
if isinstance(cos_theta_t, complex):
print('Total internal reflection')
cos_theta_t = 0
if self.polarisation == "s":
t = (2. * self.n1 * cos_theta_i) / (self.n1 * cos_theta_i + self.n2 * cos_theta_t)
else:
t = (2. * self.n1 * cos_theta_i) / (self.n1 * cos_theta_t + self.n2 * cos_theta_i)
return self.create_wave(theta=plot_theta_t, amplitude=t, material=2)
def reflect(self):
if self.n1 == self.n2:
print('Refractive indices equal - no reflection')
return None
cos_theta_i = np.cos(self.theta)
theta_r = self.theta
cos_theta_t = self.cos_snell(self.n1, self.n2, self.theta)
if isinstance(cos_theta_t, complex):
cos_theta_t = 0
plot_theta_r = -theta_r
if self.polarisation == "s":
r = (self.n1 * cos_theta_i - self.n2 * cos_theta_t) / (self.n1 * cos_theta_i + self.n2 * cos_theta_t)
else:
r = (self.n1 * cos_theta_t - self.n2 * cos_theta_i) / (self.n1 * cos_theta_t + self.n2 * cos_theta_i)
return self.create_wave(theta=plot_theta_r, amplitude=r, material=1, reverse_phase=True)
def find_critical_angle(self):
if self.n1 > self.n2:
return 0.5*np.pi - np.arcsin(self.n2/self.n1)
else:
return 0.5*np.pi + 0.0001
@staticmethod
def cos_snell(n1, n2, theta_i):
"""
Returns cosine of angle of transmission using Snell's law
Args:
n1 (float) - incident medium's refractive index
n2 (float) - transmissive medium's refractive index
theta_i (float) - angle of incidence
"""
cos_squared_theta_t = 1. - ((n1/n2) * np.sin(theta_i))**2
if cos_squared_theta_t < 0:
return cmath.sqrt(cos_squared_theta_t)
else:
return np.sqrt(cos_squared_theta_t)
def plot(self):
fig = plt.figure(figsize=(5, 5))
im_i = plt.pcolormesh(self.incident["x"], self.incident["y"], self.incident["zz"][0], vmin=-1, vmax=1)#, animated=True)
im_t = plt.pcolormesh(self.transmitted["x"], self.transmitted["y"], self.transmitted["zz"][0], vmin=-1, vmax=1)#, animated=True)
if self.frustration:
im_f = plt.pcolormesh(self.frustrated["x"], self.frustrated["y"], self.frustrated["zz"][0], vmin=-1, vmax=1)#, animated=True)
def update_plot(*args):
im_i.set_array(self.incident["zz"][self.frame][:-1, :-1].ravel())
im_t.set_array(self.transmitted["zz"][self.frame][:-1, :-1].ravel())
if self.frustration:
im_f.set_array(self.frustrated["zz"][self.frame][:-1, :-1].ravel())
self.frame += 1
if self.frame == self.num_frames:
self.frame = 0
return im_i, im_t
ani = animation.FuncAnimation(fig, update_plot, interval=50, blit=True)
plt.show()
return ani
In [ ]:
boundary = Boundary(angle=42, n1=1.5, n2=1., polarisation="p", interference=False, frustration=True)
ani = boundary.plot()